/*
 * BatchReactor.cpp
 *
 *  Created on: 3 Jun 2011
 *      Author: allan
 */

#include "BatchReactor.h"

// C++ includes
#include <cmath>
#include <functional>
#include <iostream> // TODO remove
#include <iomanip> // TODO remove
using namespace std::placeholders;

// GSL includes
#include <gsl/gsl_odeiv.h>

// FluidReact includes
#include "MultirootNewtonRaphson.h"
#include "Speciation.h"
#include "SystemOutput.h"

BatchReactor::BatchReactor() :
system(0)
{}

BatchReactor::~BatchReactor()
{}

void BatchReactor::SetSystem(System* system)
{
	this->system = system;
}

const System*const BatchReactor::GetSystem() const
{
	return system;
}

const vector<double> TotalConcentrationOfPrimarySpecies(const System* system)
{
	const unsigned numOfPrimarySpecies = system->GetNumOfPrimarySpecies();
	
	const unsigned numOfEquilibriumSpecies = system->GetNumOfEquilibriumSpecies();
	
	vector<double> Uj(numOfPrimarySpecies);
	
	for(unsigned iPrimary = 0; iPrimary < numOfPrimarySpecies; ++iPrimary)
		Uj[iPrimary] = system->GetMolesOfSpecies(system->GetIndexOfPrimarySpecies(iPrimary));
	
	for(unsigned iEquilibrium = 0; iEquilibrium < numOfEquilibriumSpecies; ++iEquilibrium)
	{
		for(unsigned iProduct = 0; iProduct < system->GetEquilibriumReaction(iEquilibrium).GetNumOfProducts(); ++iProduct)
		{
			const unsigned j = system->GetEquilibriumReaction(iEquilibrium).GetIndexOfProduct(iProduct);
			
			const double vej = system->GetEquilibriumReaction(iEquilibrium).GetStoichiometryOfProduct(iProduct);
			
			Uj[j] += vej * system->GetMolesOfSpecies(system->GetIndexOfEquilibriumSpecies(iEquilibrium));
		}
	}
	
	return Uj;
}

struct BatchReactionParameters
{
	const System* system;

	const vector<double>& ratesOfKineticSpecies;

	const vector<double>& ratesOfPrimarySpecies;
};

int ODESystemFunction(double t, const double y[], double dydt[], void* params)
{
	BatchReactionParameters p = *(BatchReactionParameters*)params;

	const unsigned numOfPrimarySpecies = p.system->GetNumOfPrimarySpecies();

	const unsigned numOfKineticSpecies = p.system->GetNumOfKineticSpecies();

	for(unsigned j = 0; j < numOfPrimarySpecies; ++j)
	{
		dydt[j] = p.ratesOfPrimarySpecies[j];
	}
	
	for(unsigned k = 0; k < numOfKineticSpecies; ++k)
	{
		dydt[k + numOfPrimarySpecies] = p.ratesOfKineticSpecies[k];
	}
	
	return GSL_SUCCESS;
}

void BatchReactor::Solve(double tEnd, double h, double hMax, string filename)
{
	const unsigned numOfPrimarySpecies = system->GetNumOfPrimarySpecies();
	
	const unsigned numOfKineticSpecies = system->GetNumOfKineticSpecies();

	//const gsl_odeiv_step_type * type = gsl_odeiv_step_rk8pd;
	const gsl_odeiv_step_type * type = gsl_odeiv_step_rk4;

	gsl_odeiv_step* s    = gsl_odeiv_step_alloc(type, numOfPrimarySpecies + numOfKineticSpecies);
	//gsl_odeiv_control* c = gsl_odeiv_control_y_new(0.0, 1.0E-10);
	gsl_odeiv_control* c = 0;gsl_odeiv_control_yp_new(1.0E-20, 1.0E-20); // limiting derivatives (todo you need to understand the adaptive algorithm better)
	gsl_odeiv_evolve* e  = gsl_odeiv_evolve_alloc(numOfPrimarySpecies + numOfKineticSpecies);

	gsl_odeiv_system sys = {ODESystemFunction, 0, numOfPrimarySpecies + numOfKineticSpecies, 0};

	double t = 0.0;
	double* y    = new double[numOfPrimarySpecies + numOfKineticSpecies];
	double* dydt = new double[numOfPrimarySpecies + numOfKineticSpecies];

	vector<double> Ck(numOfKineticSpecies);
	
	for(unsigned iKinetic = 0; iKinetic < numOfKineticSpecies; ++iKinetic)
	{
		Ck[iKinetic] = system->GetMolesOfSpecies(system->GetIndexOfKineticSpecies(iKinetic));
	}
	
	vector<double> Uj = TotalConcentrationOfPrimarySpecies(system);

	for(unsigned j = 0; j < numOfPrimarySpecies; ++j)
	{
		y[j] = Uj[j];
	}

	for(unsigned k = 0; k < numOfKineticSpecies; ++k)
	{
		y[k + numOfPrimarySpecies] = Ck[k];
	}
	
	SystemOutput sysOutput(system, "s");
	
	Speciation speciation;
	
	speciation.SetSystem(*system);
	
	for(unsigned iPrimary = 0; iPrimary < numOfPrimarySpecies; ++iPrimary)
		speciation.ImposeUj(iPrimary, Uj[iPrimary]);
	
	while(t < tEnd)
	{
		cout << fixed;
		cout << "Progress: " << t/tEnd * 100 << "%" << "\t";
		
		sysOutput.Register(t, system);
		
		vector<double> ratesOfKineticSpecies = system->GenerationRatesOfKineticSpecies();
		
		vector<double> ratesOfPrimarySpecies = system->GenerationRatesOfPrimarySpecies(ratesOfKineticSpecies);
		
		BatchReactionParameters parameters = {system, ratesOfKineticSpecies, ratesOfPrimarySpecies};

		sys.params = &parameters;
		
		int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &t, tEnd, &h, y);
		
		if(h > hMax) h = hMax;
		
		for(unsigned j = 0; j < numOfPrimarySpecies; ++j)
		{
			Uj[j] = y[j];
		}

		for(unsigned k = 0; k < numOfKineticSpecies; ++k)
		{
			Ck[k] = y[k + numOfPrimarySpecies];
		}
		
		system->SetMolesOfKineticSpecies(Ck);
		
		for(unsigned iPrimary = 0; iPrimary < numOfPrimarySpecies; ++iPrimary)
			speciation.ImposeHj(iPrimary, Uj[iPrimary]);
		
		bool converged = speciation.Solve(1.0E-6, 100);
		
		if(!converged)
		{
			sysOutput.Print(filename);
			cout << "Error: The numerical solution of the nonlinear system did not converge." << endl;
			exit(0);
		}
		
		if(status != GSL_SUCCESS)
		{
			cout << "Failure in GSL ODE solver. Check what happened, Allan." << endl;
			break;
		}
		
		cout << endl;
	}
	
	gsl_odeiv_evolve_free(e);
	gsl_odeiv_control_free(c);
	gsl_odeiv_step_free(s);

	delete [] y;
	delete [] dydt;
	
	sysOutput.Print(filename);
}
